% fig6a.m ——— 优化后的自适应 α 下 x(τ)（复现论文 Fig.6(a)，修正索引）
clear; clc; close all;

%% —— 参数设置 —— 
a     = 0.5;
b     = 0.215;
beta  = 0.2;
A     = 1.15;
omega = 0.61;

sigma1 = 0.002;       % 增益 σ₁
lambda = 1.93;        % 阈值 λ

%% —— 时间设置 —— 
dt    = 0.01;
Tmax  = 4000;         % 绘制到 τ = 1000
Nstep = round(Tmax/dt);

tau  = 0:dt:Tmax;     % 长度 Nstep+1
cosw = cos(omega * tau);

%% —— 预分配 —— 
x     = zeros(1, Nstep+1);
y     = zeros(1, Nstep+1);
z     = zeros(1, Nstep+1);
alpha = zeros(1, Nstep+1);
H     = zeros(1, Nstep+1);

% 初始条件
x(1)     = 0.01;
y(1)     = 0.02;
z(1)     = 0.02;
alpha(1) = 1.8;

%% —— 主循环（内联 RK4 + 内联能量） —— 
for k = 1:Nstep
    % 当前状态与 α
    X1 = [x(k); y(k); z(k)];
    a1 = alpha(k);

    % 四阶 RK4 的余弦值
    cw1 = cosw(k);
    cw2 = cosw(k+1);

    % —— RK4 k1 —— 
    dx1 = -a1*((X1(2)-X1(1)) - (X1(2)-X1(1))^3/3 + (a+3*b*X1(3)^2)*X1(1));
    dy1 =      (X1(2)-X1(1)) - (X1(2)-X1(1))^3/3 - beta*X1(2) + A*cw1;
    dz1 =      X1(1);
    K1  = [dx1; dy1; dz1];

    % —— RK4 k2 —— 
    X2 = X1 + 0.5*dt*K1;
    dx2 = -a1*((X2(2)-X2(1)) - (X2(2)-X2(1))^3/3 + (a+3*b*X2(3)^2)*X2(1));
    dy2 =      (X2(2)-X2(1)) - (X2(2)-X2(1))^3/3 - beta*X2(2) + A*cw2;
    dz2 =      X2(1);
    K2  = [dx2; dy2; dz2];

    % —— RK4 k3 —— 
    X3 = X1 + 0.5*dt*K2;
    % reuse cw2 as approximate cos(ω(t+dt/2))
    dx3 = -a1*((X3(2)-X3(1)) - (X3(2)-X3(1))^3/3 + (a+3*b*X3(3)^2)*X3(1));
    dy3 =      (X3(2)-X3(1)) - (X3(2)-X3(1))^3/3 - beta*X3(2) + A*cw2;
    dz3 =      X3(1);
    K3  = [dx3; dy3; dz3];

    % —— RK4 k4 —— 
    X4 = X1 + dt*K3;
    % use cw2 again as cos(ω(t+dt))
    dx4 = -a1*((X4(2)-X4(1)) - (X4(2)-X4(1))^3/3 + (a+3*b*X4(3)^2)*X4(1));
    dy4 =      (X4(2)-X4(1)) - (X4(2)-X4(1))^3/3 - beta*X4(2) + A*cw2;
    dz4 =      X4(1);
    K4  = [dx4; dy4; dz4];

    % 更新状态
    Xnew     = X1 + dt/6 * (K1 + 2*K2 + 2*K3 + K4);
    x(k+1)   = Xnew(1);
    y(k+1)   = Xnew(2);
    z(k+1)   = Xnew(3);

    % —— 内联能量计算 —— 
    H(k+1) = 0.5*a1*x(k+1)^2 + 0.5*y(k+1)^2 ...
           + 0.5*a*z(k+1)*x(k+1) + 1.5*b*z(k+1)^3*x(k+1);

    % —— 自适应 α 更新 —— 
    if H(k+1) - H(k) > lambda
        alpha(k+1) = a1 + sigma1 * a1 * dt;
    else
        alpha(k+1) = a1;
    end
end

%% —— 绘图 —— 
figure;
plot(tau, x, 'm', 'LineWidth',1.2);
xlabel('\tau','FontSize',12);
ylabel('x','FontSize',12);
title('(a)','FontSize',14,'FontWeight','normal');
xlim([0 Tmax]); 
grid on;

